Matlab
code 2.3: Matlab file “Figure 2-7.m”
%--------------------------------------------------------------------
% This code can be used to generate Figure 2.7
%--------------------------------------------------------------------
clear all;
close all;
%% the theoretical curve of precession target
micro-Doppler characteristic
c = 3e8;
j = sqrt(-1);
fc = 10e9; % carrier frequency of transmitted
signal
v = 0; % translational velocity of target
cord = 1000*[300 100 500]; % coordinates of
local coordinate system's origin in the radar coordinate system
colo = [0 0 0]; % coordinates of radar in the
radar coordinate system
R0 = cord-colo;
n = R0/sqrt(sum(R0.^2)); % unit vector of
Line-of-Sight(LOS)
p1 = [0 0 0.5]; % scatterer on the apex of cone
p2 = [0.5 0 -0.5]; % scatterer on the bottom of
cone
tgt = [p1;p2];
ae = pi*[1/3 1/4 1/5]; % Euler angle
ws = [19.2382 -11.1072 22.2144]; % angular
velocity of spinning in the local coordinate system
omegas = 10*pi; % angular frequency of spinning
% omegas = sqrt(sum(ws.^2));
Ts = 0.2; % spinning period
wc = [7.6756 -2.3930 9.6577]; % angular
velocity of coning in the local coordinate system
omegac = 4*pi; % angular frequency of coning
% omegac = sqrt(sum(wc.^2));
Tc = 0.5; % coning period
t = 2; % radar illumimated time
ri1 = [cos(ae(3)) sin(ae(3)) 0;-sin(ae(3))
cos(ae(3)) 0;0 0 1];
ri2 = [1 0 0;0 cos(ae(2)) sin(ae(2));0
-sin(ae(2)) cos(ae(2))];
ri3 = [cos(ae(1)) sin(ae(1)) 0;-sin(ae(1))
cos(ae(1)) 0;0 0 1];
ri = ri1*ri2*ri3; % initial rotation matrix
ws = ri*ws'; % angular velocity of spinning in
the reference coordinate system
% omega = sqrt(sum(w.^2));
wse = ws/omegas; % unit vector of spinning
wsr = [0 -wse(3) wse(2);wse(3) 0
-wse(1);-wse(2) wse(1) 0]; % skew symmetric matrix of spinning
wc = ri*wc'; % angular velocity of coning in
the reference coordinate system
wce = wc/omegac; % unit vector of coning
wcr = [0 -wce(3) wce(2);wce(3) 0
-wce(1);-wce(2) wce(1) 0]; % skew symmetric matrix of coning
r0a = tgt(1,:)-colo;
r0b = tgt(2,:)-colo;
r0ar = ri*r0a'; % spinning scatterer in the
reference coordinate system
r0br = ri*r0b'; % coning scatterer in the
reference coordinate system
a1 = (2*fc/c)*r0br'*(omegac*wcr^2+omegac*wcr^2*wsr^2)'*n';
a2 =
(2*fc/c)*r0br'*(omegac*wcr+omegac*wcr*wsr^2)'*n';
a3 =
(2*fc/c)*r0br'*(omegac*wcr^2*wsr+omegas*wcr*wsr^2)'*n';
a4 =
(2*fc/c)*r0br'*(omegas*wcr*wsr-omegac*wcr^2*wsr^2)'*n';
a5 =
(2*fc/c)*r0br'*(omegac*wcr*wsr-omegas*wcr^2*wsr^2)'*n';
a6 =
(2*fc/c)*r0br'*(-omegac*wcr*wsr^2-omegas*wcr^2*wsr)'*n';
a7 =
(2*fc/c)*r0br'*(omegas*wsr^2+omegas*wcr^2*wsr^2)'*n';
a8 =
(2*fc/c)*r0br'*(omegas*wsr+omegas*wcr^2*wsr)'*n';
prf = 4000; % pulse repetition frequency
pri = 1/prf; % pulse repetition interval
dt = 0:pri:t-pri; % time sampling interval
m = length(dt);
fmd = zeros(2,m); % theoretical micro-Doppler
frequency
for i = 1:m
fmd(1,i) =
a1*sin(omegac*dt(i))+a2*cos(omegac*dt(i))+a3*sin(omegas*dt(i))*sin(omegac*dt(i))...
+a4*sin(omegac*dt(i))*cos(omegas*dt(i))+a5*cos(omegac*dt(i))*sin(omegas*dt(i))...
+a6*cos(omegac*dt(i))*cos(omegas*dt(i))+a7*sin(omegas*dt(i))+a8*cos(omegas*dt(i));
fmd(2,i) =
(2*fc*omegac/c).*((wcr^2.*sin(omegac*dt(i))+wcr.*cos(omegac*dt(i)))*ri*r0ar)'*n';
end
figure(1)
plot(dt,fmd)
xlabel('Time (s)')
ylabel('Frequency (Hz)')
axis([0,2,-1500,1500])
%% the time-frequency analysis result of
precession target micro-Doppler characteristic
r = zeros(length(tgt(:,1)),length(dt)); %
distance between the scatterers and radar
for i = 1:m
for n
= 1:length(tgt(:,1))
if n == 1
r(n,i) =
sqrt((R0'+v*dt(i)+((eye(3)+wcr*sin(omegac*dt(i))+wcr^2*(1-cos(omegac*dt(i))))*r0ar))'...
*(R0'+v*dt(i)+((eye(3)+wcr*sin(omegac*dt(i))+wcr^2*(1-cos(omegac*dt(i))))*r0ar)));
else
r(n,i) =
sqrt((R0'+v*dt(i)+(eye(3)+wcr*sin(omegac*dt(i))+wcr^2*(1-cos(omegac*dt(i))))...
*(eye(3)+wsr*sin(omegas*dt(i))+wsr^2*(1-cos(omegas*dt(i))))*r0br)'...
*(R0'+v*dt(i)+(eye(3)+wcr*sin(omegac*dt(i))+wcr^2*(1-cos(omegac*dt(i))))...
*(eye(3)+wsr*sin(omegas*dt(i))+wsr^2*(1-cos(omegas*dt(i))))*r0br));
end
end
end
s = zeros(length(dt),length(tgt(:,1))); % echo signal
matrix
for i = 1:length(tgt(:,1))
s(:,i) = exp(j*2*pi*fc*2*r(i,:)'/c);
end
st = sum(s,2);
N = 200; % number of Gabor coefficients in time
Q = 50; % degree of oversampling
tfr = tfrgabor(st.',N,Q); % Gabor transform
tfr_r = fftshift(tfr,1);
figure(2)
contour(linspace(min(dt),max(dt),length(tfr_r(1,:))),linspace(-prf/2,prf/2-1,length(tfr_r(:,1))),tfr_r,70)
xlabel('Time (s)')
ylabel('Frequency (Hz)')
axis([0,2,-1500,1500])